Averaging For Solitons With Nonlinearity Management 
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We develop an averaging method for solitons of the nonlinear Schrodinger equation with pe- 
riodically varying nonlinearity coefficient. This method is used to effectively describe solitons in 
Bose-Einstein condensates, in the context of the recently proposed and experimentally realizable 
technique of Feshbach resonance management. Using the derived local averaged equation, we study 
matter-wave bright and dark solitons and demonstrate a very good agreement between solutions of 
the averaged and full equations. 



Introduction. Dispersive nonlinear wave equations are 
appropriate mathematical models for various nonlinear 
phenomena in fluid mechanics, optics, plasmas and con- 
densed matter physics. The prototypical equation that 
generically emerges in the description of envelope waves 
is the nonlinear Schrodinger (NLS) equation [1-3] of the 
form: 



iu t = —DAu 



T\u\ 2 u- 



V(x)u. 



(1) 



Here u(x, t) is a complex envelope field, V(x) is an ex- 
ternal potential, A is the Laplacian operator in multi 
dimensions, and D and V are the coefficients of the dis- 
persive and nonlinear terms respectively. 

In a number of physical applications, the coefficients 
D and T exhibit temporal periodic variations. When 
D = D(t), the NLS equation (1) describes the dispersion 
management (DM) scheme in fiber optics, which is based 
on periodic alternation of fibers with opposite signs of 
the group-velocity dispersion. The DM scheme supports 
robust breathing solitons [4], which are well described 
through the averaging method by the integral NLS equa- 
tion [5] . Extensions of the averaging method were devel- 
oped for strong management with large variations of the 
dispersion coefficient [6] and for weak management with 
small variations of the dispersion coefficient [7]. 

When r = L(i), the NLS equation (1) has applications 
in optics for transverse beam propagation in layered op- 
tical media [8] , as well as in atomic physics for the Fesh- 
bach resonance [9] of the scattering length of inter-atomic 
interactions in Bose-Einstein condensates (BECs). The 
periodic variation of the scattering length by means of an 
external magnetic field provides an experimentally real- 
izable protocol for the generation of robust matter-wave 
breathers [10], and for their persistence against collapse 
type phenomena in higher dimensions [11,12]. Solitary 
waves have become a focal point in studies of BEC both 
theoretically and experimentally [13,14] due to their co- 
herence properties. Hence, nonlinearity management us- 
ing Feshbach resonance promises to provide a viable al- 
ternative for the generation of coherent nonlinear wave 
structures. 



Given the importance in nonlinear optics and con- 
densed matter physics, of applications of the NLS equa- 
tion (1) with periodically varying nonlinearity coefficient, 
we extend the averaging method of [5,6] to solitons with 
strong nonlinearity management, when the periodic vari- 
ations of the nonlinearity coefficient are large in ampli- 
tude. Comparing with earlier works, we note that the 
averaged equation for strong dispersion management in 
[5,6] is nonlocal, whereas our main averaged equation 
(see Eq. (10)) for strong nonlinearity management is 
local. Furthermore, our averaging method is more gen- 
eral than the asymptotic expansion method, exploited for 
weak dispersion management in [7] and for weak nonlin- 
earity management in [11]. Since the averaged equation 
obtained herein is simple, we compute numerically soli- 
tary waves of the averaged equation and compare with 
those of the full problem, showing the excellent agree- 
ment between the two. 

We emphasize that the main contribution of this work 
is two-fold. From a mathematical point of view, it is the 
derivation of a novel averaged equation that describes dy- 
namics of solitary waves under nonlinearity management. 
From a physical point of view, the main result is the 
computation of the parameter domains, where nonlinear 
waves exist in the BEC and nonlinear optics. We hope 
that highlighting the relevant analogies and differences 
may also stimulate additional cross-fertilization between 
these sub-disciplines and their respective mathematical 
techniques. Furthermore, our work poses the interesting 
problem of understanding what happens in the parameter 
domains, where solitary waves do not exist. These fun- 
damental problems are of interest not only to theorists 
and experimentalists in atomic and optical physics but 
also, more generally, to researchers in nonlinear and wave 
physics, where periodic temporal variations and their av- 
eraging methods are studied (see e.g., [15]). 

Derivation of the averaged equation. We start with the 
NLS equation (1) with D = 1 and T = T(t). The exter- 
nal potential V(x) is left arbitrary but we keep in mind 
that the magnetic and laser trappings relevant to BEC 
applications impose parabolic and periodic potentials re- 
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spectively. Also, we will restrict ourselves to one spatial 
dimension, but generalization of the method to multi- 
dimensions is straightforward. The resulting equation 
(also referred to as the Gross-Pitaevskii (GP) equation 
[3]) describes the "cigar-shaped" BECs and reads: 

iu t = -u xx + T(t)\u\ 2 u + V{x)u, (2) 

where the nonlinearity coefficient (proportional to the 
scattering length in the BECs) T(t+e) = T(t) is a smooth, 
sign-indefinite, periodic function of period e. We assume 
that the period e of the nonlinearity management is small 
compared to the characteristic propagation time of non- 
linear waves, while the nonlinearity variations are large 
in amplitude. In this case, we decompose T(i) into the 
mean- value part 70 and a large fast- varying part 7, ac- 
cording to the representation: 

r(*)= 7 o + ^7(r), r=- e , (3) 

where j(t + 1) = j(t) and j(r)dr = 0. Using 

u(x,t) = v(x,r)exp ^-i J 7(t')|v| 2 (x, t')c?t'^ , (4) 

we remove the large fast variations of the nonlinearity co- 
efficient and bring the GP equation (2) to an equivalent 
form: 

ie~ 1 v T = -v xx + -f \v\ 2 v + V(x)v 

pT pt 

+2iv x j(r')\v\ 2 x (x,T')dT' + iv j(r')\v\l x (x,T')dT' 
Jo Jo 

+v^J\(r')\v\U Xl r')dT'J . (5) 

In the averaging method (see [6] for details), we decom- 
pose solutions of the problem with variable coefficients 
(5) into a slowly varying mean part w(x,t) and a small, 
fast- varying part Vi(x, t): 

v(x, t) = w(x, t) + evi(x, r; w(x, t)), t = er. (6) 

The varying part v\(x,t;w) is a periodic function of r 
with unit period. To leading order, this condition is sat- 
isfied if w(x,t) satisfies the averaged equation: 

iw t — - w xx + 7o|w| 2 w + V(x)w 

+ 2iv 1 w x \w\l + iv!w\w\l x + v 2 w {\w\ 2 x ) 2 , (7) 

where v\ — j^v(r)dT, v 2 = J^v 2 {r)dT, and v(t) = 
Jo 7( r ')^ T '- The averaging method is simplified with the 
gauge transformation, 

w(x,t) — ip(x,t)exp (ii/i\tp\ 2 (x,t)) , (8) 

which reduces (7) to the form: 



iipt - viip\ip\t = -^xx+ 7o|V'| 2 ^ + V(x)tp 

+ (9) 

where /i = v 2 — v\. Using the balance equation i|V'| 2 = 
('ipxip — 'ipipx) x , which follows from (9), we rewrite the 
averaged equation in the final form: 

#t = - rpxx + loH\ 2 4> + V{x)il> + MV> (IV'Ix) 2 

+ iviip{i}i) x -i) x if) x . (10) 

The averaged equation (10) is the main result of this 
Letter. It is seen to be equivalent to the integral aver- 
aged equation derived for strong dispersion management 
in fiber optics [5,6], but it is a local evolution equation. A 
similar local equation was also derived for weak disper- 
sion management in fiber optics [7], when the last two 
terms of (10) are small compared to the leading-order 
NLS equation. We emphasize that our main averaged 
equation (10) is derived for strong nonlinearity manage- 
ment and it captures all terms in the same, leading order 
of the averaging method. 

Solitons in BECs. The simplest standing waves of the 
averaged equation (10) are obtained through the stan- 
dard ansatz [1]: 

4>(x,t) = 4>{x)e i » t , (11) 

where 4>{x) solves the second-order differential equation: 

-4>" + w<f> + V{x)(j) + 7o 3 + W) 2 3 = 0. (12) 

As a typical example of a smooth periodic variation of 
the scattering length [11,12], we use the sinusoidal func- 
tion T(t) = 70 +71 sin(27rf), in which case fj, — 7 2 /(87r 2 ). 
We also set e = 1 and choose \u>\ € [0.1,0.5] to ensure 
validity of the averaged equation (12), when e <C 27r/|u;|. 
We also use the parabolic potential V(x) for the mag- 
netic trapping of the BEC, V(x) — \Vl 2 x 2 , where fi 2 e 
[0.02,0.4]. 

To estimate actual physical quantities corresponding 
to the above values of the normalized parameters, we 
first note that the cases 70 < (70 > 0) are relevant to 
an attractive (repulsive) condensate, such as 7 Li ( 85 Rb), 
characterized by a negative (positive) scattering length 
a = — lnm (a = 0.8nm), in a magnetic field B w 650 
G (B « 159 G). These values of the scattering lengths 
set the units in the parameters 70 and 71, which may 
take different values as long as the magnetic field B is 
varied [9]. On the other hand, the number of atoms 
N in the two condensates is taken to be as follows: 
N = 1 x 10 4 (N = 2 x 10 5 ) for fl 2 = 0.4 (fl 2 = 0.02) for 
the 7 Li condensate and N = 4 x 10 3 (N = 7.5 x 10 4 ) for 
O 2 = 0.4 (n 2 = 0.02) for the 85 Rb condensate. Since we 
deal with cigar-shaped BECs, we may assume that the 
external magnetic trap is highly anisotropic, character- 
ized by the confining frequencies u\\ — 2tt x 3.6Hz and 
uj± = 2-7T x 360Hz in the axial and transverse directions 
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respectively. In such a case, the time and space units in 
the results that will be presented below, are 44.2ms and 
2/j,m (for 7 Li) or 44.2ms and 0.6^m (for 85 Rb). 

Numerical Results. Using Eq. (12), we can now ob- 
tain the solution 4>{x) for a given set of parameters 
(70, 71, 0, u), by means of the Newton method. We also 
perform parameter continuations, to follow the solution 
branches as the parameters vary. 

Fig. 1 shows two solutions of the averaged equation 
(12) with 70 < (the attractive BEC with negative scat- 
tering length), 71 = 0.5, il 2 = 0.4, and u = 0.5. The so- 
lution on the left panel is the bright soliton, which has the 
form 4>{x) = (2cj/7o) 1//2 sech(a-> 1 / 2 a;) when 71 = ft = 0. 
The solution on the right panel is the so-called twisted 
soliton, which corresponds to a concatenation of two sep- 
arated bright solitons of opposite parity (see e.g., [16]). 
The twisted soliton does not exist when 70 < and 
CI = 0. Higher-order solutions with multiple nodes (ze- 
ros) may also exist in the averaged equation (12) with 
70 < and fl =/= 0, in some parameter domains. 

Fig. 2 shows two solutions of the averaged equation 
(12) with 70 > (the repulsive BEC with positive scat- 
tering length), 71 = 0.5, and uj = —0.5. In the case of 
70 > 0, the localized solutions of Eq. (12) bifurcate from 
linear modes trapped by the parabolic potential V(x), 
such that an infinite number of solitons with multiple 
nodes (zeros) exists for larger negative values of uj. The 
solution on the left panel for £1 2 = 0.4 is the ground 
state, often approximated by the Thomas-Fermi cloud 
[10]. The solution on the right panel for CI 2 = 0.02 is 
the (embedded in the Thomas-Fermi cloud) dark soli- 
ton which, in the case of 71 = £1 = 0, has the form 
cj){x) = (|w|/7o) 1/2 tanh((|w|/2) 1 / 2 x) when 71 = ft = 0. 
The dark soliton is the only localized solution of Eq. (12) 
with 70 > and O = 0. Notice that the regular dark 
soliton asymptotes to a non-vanishing amplitude when 
£7 = 0, while it asymptotes to in the presence of the 
magnetic trap. 

Fig. 3 shows two parameter (70, 71) continuation of the 
dark (top panel for f2 2 = 0.02) and bright (bottom panel 
for fi 2 = 0.4) soliton solutions of the averaged equation 
(12) with \uj\ = 0.5. The branch of dark solitons exists 
above the bifurcation curve on the top panel, whereas 
the branch of bright solitons exists below the curve on 
the bottom panel. The two curves pass through the ori- 
gin 70 = 71 = 0. We note that the domain of existence 
of dark and bright solitons shrinks for increasing values 
of 71. 

Finally, we examine how well the averaged equation 
(12) approximates bright and dark solitons of the GP 
equation (2). In our numerical simulations of Eq. (2), 
we initialize the wavefunction, using the spatial profile 
obtained from (12), and subsequently observe whether 
the temporal evolution of Eq. (2) preserves the average 
profile of Eq. (12). 

Fig. 4 shows the temporal evolution of the bright soli- 



ton with 70 = —0.5, 71 = 1, O 2 = 0.4, and to = 0.5. 
The periodic variations of the nonlinearity coefficient T(t) 
in Eq. (2) lead to complicated oscillations of the solu- 
tion's maximum. While the solution oscillates between 
single-humped and double-humped solitons (a scenario 
that bears analogies to the observations of [10]), the av- 
erage of the two extreme solitons (at the maximum and 
the minimum amplitudes) is practically indistinguishable 
from the profile <j>(x) of Eq. (12). 

Fig. 5 shows the temporal evolution of the dark soli- 
ton with 70 = 0.5, 71 = 1, fl 2 = 0.02, and u = -0.5. 
We notice that the center of the dark soliton remains 
at the origin x — 0, without any oscillations. Only the 
maxima of |u| 2 (x, t) display periodic oscillations of small 
amplitude. The average of the extreme solitons is again 
essentially identical to the profile 4>(x) of the averaged 
equation (12). 

Conclusion. We have derived and studied the aver- 
aged equation (10) for the NLS (GP) equation (2) with 
periodic modulation of the nonlinearity coefficient. Our 
results are of broad interest to diverse areas of atomic 
and optical physics, as well as of nonlinear and, also, 
mathematical physics. We have identified numerically 
several branches of solitary waves of the averaged sta- 
tionary equation (12). We have also compared solutions 
of the averaged and full equations, obtaining a very good 
agreement. It is of particular and immediate interest to 
examine these predictions experimentally, as well as to 
identify what happens in parametric regimes where the 
present theory is no longer able to identify such waves. 
Furthermore, the averaged equation (10) can be useful for 
future studies of BECs, e.g., in optical lattice potentials 
or in higher dimensions. 

This work was supported in part by NSERC grant # 
RGP238931 (DEP), and by a UMass FRG, NSF-DMS- 
0204585 and the Eppley Foundation for Research (PGK). 
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FIG. 1. Bright solitons (left) and twisted solitons (right) obtained from Eq. (12) with 70 < 0, 71 = 0.5, Q. 2 = 0.4, and 
ui = 0.5. The top subplot shows the solution maximum for different values of 70. The bottom subplot shows the potential 
(dashed line) and the solutions: at the left panel the solution is shown for 70 = —0.8 (solid), —0.6 (dash-dotted), —0.4 (dotted) 
—0.2 (circles) and —0.02 (stars) and at the right panel for 70 = —0.8 (solid), —0.6 (dash-dotted) and —0.37 where the branch 
ceases to exist (circles). 




FIG. 2. Similar to the first figure but for Thomas-Fermi clouds (left) and dark solitons (right) obtained from Eq. (12) 
with 70 > 0, 71 = 0.5, u — 0.5, while Q 2 = 0.4 (left) and Q. 2 = 0.02 (right). The top subplot shows the solution maximum 
for different values of 70. The bottom subplot shows the potential (dashed line) and the solutions for 70 = 0.8 (solid) , 0.6 
(dash-dotted), 0.4 (dotted), 0.2 (circles) and 0.01 (stars). Notice that in the case of the dark soliton in the right subplot, the 
solution profile is altered very slightly between the cases 70 = 0.8 and 70 = 0.2, which are practically indistinguishable. 
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FIG. 3. Domain of existence of dark (top panel for Q 2 = 0.02) and bright (bottom panel for Q, 2 = 0.4) solitons of Eq. (12) 
with \uj\ = 0.5. The solutions exist above and below, respectively, the corresponding curves of the (70,71) plane. 
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FIG. 4. Temporal evolution of the bright soliton with 70 = —0.5, 71 = 1, Q. 2 — 0.4, and uj = 0.5. The top panel shows 
evolution of the (spatial) maximum of |u| 2 (a;, t) as a function of time. The middle panel shows the solution |u| 2 (a;, t) at t = 137.3 
(solid) and t = 138 (dashed), their average (circles), and the initial configuration (dash-dotted). The latter practically coincides 
with the average. The bottom panel shows a contour plot of |u| 2 (a;,t) in (x,t). 
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